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ABSTRACT 

We investigate the evolution of protoplanets with different masses embedded in an accretion disk, 
via global fully three-dimensional hydrodynamical simulations. We consider a range of planetary 
masses extending from one and a half Earth's masses up to one Jupiter's mass, and we take into 
account physically realistic gravitational potentials of forming planets. In order to calculate accu- 
rately the gravitational torques exerted by disk material and to investigate the accretion process 
onto the planet, the flow dynamics has to be thoroughly resolved on long as well as short length 
scales. We achieve this strict resolution requirement by applying a nested-grid refinement technique 
which allows to greatly enhance the local resolution. Our results from altogether 51 simulations 
show that for large planetary masses, approximately above a tenth of the Jupiter's mass, migration 
rates are relatively constant, as expected in type II migration regime and in good agreement with 
previous two-dimensional calculations. In a range between seven and fifteen Earth's masses, we find 
a dependency of the migration speed on the planetary mass that yields time scales considerably 
longer than those predicted by linear analytical theories. This property may be important in de- 
termining the overall orbital evolution of protoplanets. The growth time scale is minimum around 
twenty Earth-masses, but it rapidly increases for both smaller and larger mass values. Significant 
differences between two- and three-dimensional calculations are found in particular for objects with 
masses smaller than ten Earth-masses. We also derive an analytical approximation for the numeri- 
cally computed mass growth rates. 

Subject headings: accretion, accretion disks — hydrodynamics — methods: numerical — planetary systems: 
formation 

Introduction 



Over the past few years the interest of the astronomical 
community in the enigma of planet formation and evolu- 
tion has been rising as the number of observed extrasolar 
planets has continued to increase. Today about 100 ex- 
trasolar planets are known, mostly orbiting around main- 
sequence stars of solar type. An always up-to-date status 
of extrasolar planet detections can be found at the Extra- 



solar Planet Encyclopedia , maintained by J. Schneider, 
or at the California & Carnegie Planet Search^. 

Several of the orbital and physical key properties of 
planets (e.g., location, eccentricity, rotation rate, mass) 
are believed to originate from the early phases of planet 
formation, when the protoplanet is still embedded in the 
surrounding protostellar disk from which it generated. In 
particular, the small semi-major axis of several (51 Peg- 
type) planets is usually interpreted as a migration process 
produced by gravitational torques of the disk material act- 



^To appear in The Astrophysical Journal (v586 nl 
March 20, 2003 issue). Also available as ApJ preprint doi: 
10.1086/367555. 



^http: / /www. obspm.fr/planets. 
^http:/ /exoplanets.org/. 
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ing on the protoplanets. In order to properly take these 
effects into account, one has to consider the joint evolution 
of the circumstellar disk and the embedded planet. 

Following the first and mostly linear analytical stud- 
ies (Goldreich & Tremaine 1980; Papaloizou & Lin 1984; 
Ward 1986), fully non-linear numerical sirimlations have 
been performed lately (Kley 1999; Lubow, Seibert, & 
Artymowicz 1999; Nelson et al. 2000; Papaloizou, Nel- 
son, & Masset 2001; Kley, D'Angelo, & Henning 2001; 
D'Angelo, Henning, & Kley 2002; Tanigawa & Watanabe 
2002) in order to achieve a deeper insight into the phys- 
ical processes governing the interactions between a pro- 
tostellar disk and an embedded protoplanet. Among the 
main discoveries we can cite: t) the creation of spiral den- 
sity perturbations in the disk; ii) the formation of a deep 
annular gap along the orbit of Jupiter-mass planets; in) 
an inward migration resulting from the net gravitational 
torques caused by inner and outer density wave perturba- 
tions; iv) the continuation of mass accretion through the 
gap- 
So far most of the computations have investigated 
mainly the effects due to large-scale interactions, much 
larger than the size of the Roche lobe of the forming 
planet. Furthermore, very little is known about the non- 
linear effects that such interactions have when low-mass 
planets are involved. The reason for this has been primar- 
ily the lack of appropriate rmmerical tools. 

In the majority of the previous studies, the disk is 
modeled as a two-dimensional (r-ip) system, by using 
vertically- averaged quantities. Two main arguments lie 
behind this choice. First, on a physical basis, the validity 
of a two-dimensional (2D) description is consistent be- 
cause the Hill radius of a massive object is larger or com- 
parable to the disk semi-thickness. In fact, this basically 
means that the sphere of gravitational influence of the 
embedded body, i.e., the Hill sphere'*, contains the whole 
vertical extent of the disk. But this usually implies that 
the planet must have a mass on the order of one Jupiter- 
mass. Second, a less massive planet has a weaker impact 
on the disk, requiring a higher resolution to compute prop- 
erly and highlight its effects. Such requirement typically 
rules out a full three-dimensional treatment. Although 
there is still a lot of information to be gained by per- 
forming 2D simulations (e.g., to study radiative effects or 
multi-planet systems) in particular in the case of large and 
medium mass planets, provided that the local resolution 
around the planet is accurate enough, three-dimensional 
(3D) effects become more and more important as the mass 
of the simulated planet is reduced. 

Yet, in many instances, the use of the two-dimensional 

approximation is merely dictated by the computational 
costs of 3D calculations which are generally not afford- 



^Thc Hill sphere is a measure of the volume of tfie Roehe lobe. 
In a purely restrieted three-body problem, pareels of matter in- 
side the Roche lobe are gravitationally bound to the secondary. 
Hence, they are confined to that volume of space. 



able. Depending on the resolution one is interested in, 3D 
runs still take more than an order of magnitude of CPU 
time with respect to 2D runs. As a proof of the severe 
limitations posed by fully three-dimensional calculations, 
only very few papers have been published on this issue. 

Miyoslii et al. (1999) made a comparison study of 2D 
and 3D disks within the framework of the shearing sheet 
model. Hence, rather than considering the whole disk, 
they were restricted to local sirrmlatioris. Global three- 
dimensional simulations were performed for the first time 
by Kley, D'Angelo, & Henning (2001, hereafter Paper I) 
who also measured the gas accretion rate onto the planet 
and the gravitational torques which cause the planet to al- 
ter its orbit. They found that for planetary masses below 
one half of Jupiter's mass, the outcomes of 3D calcula- 
tions start to differ from those of 2D ones. They also 
pointed out that to obtain more reliable results, the flow 
within the Roche lobe needs to be accurately resolved. 
This was done recently, for an inflnitesimally thin disk, 
by D'Angelo, Henning, & Kley (2002, hereafter Paper II) 
who introduced, for the first time in this context, a nested- 
grid technique in order to model in detail a variety of plan- 
etary masses, spanning from Earth's to Jupiter's. The 
authors proved such approach to adapt comfortably to 
these computations because global-scale structures as well 
as small- and very small-scale features of the flow can be 
captured simultaneously. They demonstrated that disks 
form around high- and low-mass planets and that circum- 
planetary material can exert very strong torques on the 
planet, usually slowing down their inward drifting motion. 

In the present paper we intend to combine the fully 
three-dimensional and global treatment of disk-planet in- 
teractions with a nested-grid refinement technique in or- 
der to carry out an extensive study on migration, accre- 
tion, and flow features around large- and small-size pro- 
toplanets. Thus, the paper comes as an extension to Pa- 
per I and Paper II. In addition, here we abandon the 
standard approach of treating the planet as a point-mass 
but rather assume that it has an extended structure. 

The outline of the paper is as follows. Section 2 deals 
with those aspects of the physical description that we 
adopt and which were not already specified in Paper I. 
We explain how we approximate the protoplanet 's struc- 
ture by using different solutions for the gravitational po- 
tential. Section 3 presents a brief overview about the nu- 
merical procedures employed in this work and describes 
the technical details of the models. As for the implemen- 
tation of the nested grids in three dimensions, for brevity 
we mainly refer to the two-dimensional strategy traced in 
Paper II. The various results of our simulations are ad- 
dressed in § 4. Fluid circulation, gravitational torques, or- 
bital migration, mass accretion rates, and how all of them 
depend on the perturbcr mass are examined. A compar- 
ison between 2D and 3D models is also carried out, to- 
gether with am analysis of some numerical effects. In § 5 
two issues related to 2D and 3D geometry effects are dis- 
cussed in more detail. We finally present our conclusions 
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in § 6. 

2. Physical Description 

The nature of most astrophysical objects is such that 
their behavior can be approximated to that of fluids. This 
is indeed the case for circumstellar disks, hence we can 
rely on the hydrodynamic formalism to describe them. 
The equations of motion that govern the evolution of a 
disk in a spherical polar coordinate system {O; -R, 0, (p} 
are presented in Paper. I and therefore, for the sake of 
brevity, we refer the interested reader to it. We assume 
the disk to be a viscous medium and include the viscos- 
ity terms explicitly by employing a complete stress tensor 
for Newtonian fluids (see e.g., Mihalas & MihaJas 1999, 
Chapter 3). 

The set of equations for the hydrodynamic variables 
{p,UR,U0,Uip) is written with respect to a reference frame 
rotating at a constant rate Q, around the polar axis 6 = 0, 
and whose origin O resides in the center of mass of the 
star-planet system. The planet is maintained on a fixed 
circular orbit, lying in the midplane of the disk (d = it/2). 
If we let Q coincide with the angular velocity of the planet 
Op, the planet docs not move within the reference frame. 
The assumption that a single protoplanet, not heavier 
than Jupiter, moves on a circular orbit is reasonable be- 
cause the global effect of the resonances, arising from disk- 
planet interactions, in most of the cases favors an eccen- 
tricity damping (Papaloizou, Nelson, & Masset 2001; Ag- 
nor & Ward 2002). 

Disk material evolves under the combined gravitational 
action of a star and a massive body. In fact, as long as the 
inner parts of low-mass protostellar disks are concerned, 
self-gravity can be neglected. Indicating with R^^ the ra^- 
dius vector pointing to the position of the star, the grav- 
itational potential $ of the whole system is represented 
by 



(1) 



where M^f is the stellar mass. In equation (1), the func- 
tion $p identifies the perturbing potential of the planet, 
which we leave unspecified for the moment. 

Since the energy equation is not considered in the 
present work, we join an established trend (e.g., Kley 
1999; Lubow, Seibert, & Artymowicz 1999; Miyoshi et 
al. 1999; Nelson et al. 2000; Papaloizou, Nelson, & Mas- 
set 2001; Tanaka, Takeuchi, & Ward 2002; Masset 2002; 
Tanigawa & Watanabe 2002) and use a locally isothermal 
equation of state as closure of the hydrodynamic equations 
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where the sound speed Cs equals the Kcplcrian velocity vk 
times the disk aspect ratio h = H/{R sin 8). The length H 
is the pressure scale-height of the disk, that also represent 
its semi-thickness. As the ratio h is assumed to be con- 
stant, the disk is azimuthally and vertically isothermal, 
whereas radially T oc h^Mi^/{R sin 9). This simplified 



approach permits to circumvent the difficulties posed by 
the solution of a complete energy equation which nobody 
has tackled yet. In fact this kind of computations would 
require a length of time which is presently not affordable. 
As reference, even without including energetic aspects, the 
CPU-time consumed by our three-dimensional global sim- 
ulations is already between ten and twenty times as long 
as that spent by two-dimensional ones. An investigation 
into the effects that may arise in two-dimensional disks 
when an energy equation is also taken into account, will 
be presented in a forthcoming paper. 

However, an important issue to improve the physical 
description of the system in the vicinity of the protoplanet 
is to adopt an appropriate equation of state which can ac- 
count for the protoplanetary envelope. Yet, in our case 
this would imply that either p or p should be specified in 
some volume around the planet. In order to avoid this, we 
choose to constrain the local structure by means of suit- 
able analytic expressions for $p . We assume that the pro- 
toplanet has a measurable size, i.e., it can be resolved by 
the employed computational mesh. Within the planetary 
volume, we approximately take into account the effects 
due to self-gravity by imposing a certain gravitational 
field. Since we aim at covering various possible scenar- 
ios, we utilize four different forms of planet gravitational 
potential, each representing a protoplanet with different 
characteristics. It is worthwhile to point out that with this 
choice none of the hydrodynamic variables {p,UR,U0,u^) 
is prescribed in any case. They simply evolve in a partic- 
ular gravitational field. Therefore, planetary material is 
allowed to interact with the surrounding environment so 
that their mutual evolutions are still connected. 

2.1. Planet Gravitational Potential 

With no exception, both numerical and analytical work 
that have so far investigated the interactions between mas- 
sive bodies and protostellar disks have made the point- 
mass assumption, i.e., the protoplanet has a finite mass 
Mp but no physical size, as was done in Paper I and 
II. This property is expressed through the gravitational 
potential 



LR — Rr, 



(3) 



where Rp is the radius vector indicating the position of 
the planet. 

Because of the singularity at i? = Rp, a parameter e is 
introduced in order to smooth the function over a certain 
region. If we denote S = R — Rp the position vector 
relative to the planet, the smoothed point-mass potential 
can be written in the following form 



1 + 



(4) 



A physical meaning of the smoothing length can be de- 
duced from equation (4). The potential <I>p enters the 
Navier-Stokes equations through its derivatives, which can 
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be reduced to d^p/dS because of the spherical symmetry 
of the gravitational field. Restricting to distances S < e, 
a binomial expansion of that derivative yields 



dS 



GMp 



GMp 



1 - 



(5) 



As we will see in § 2.1.1, equation (5) can be inter- 
preted as the sign-reversed gravitational force, per unit 
mass, exerted by a spherically homogeneous medium of 
radius e. Thus, the smoothing may act as an indicator of 
the size of the planet. 

Equation (4) may be appropriate to describe the solid 
core of a protoplanet which does not possess any signifi- 
cant envelope. Indicating with q the planet-to-star mass 
ratio Mp/Mir, in these computations we apply the po- 
tential 'I'p''' in models with g > 3 x 10^", i.e., ten times 
the Earth's mass (Me) when Al* = 1 Mq. In all cases 
the smoothing parameter e is set to 10% of the Hill ra- 
dius Rh^. Along with equation (4), we introduce three 
alternative expressions for $p, namely the potential of a 
homogeneous sphere (<I>p^); that describing a fully radia- 
tive and static envelope ($p^); and finally that proper for 
a fully convective and static envelope ($p^). A compar- 
ative example of their behavior is reported in Figure 1. 

Among the last three solutions, only the gravitational 
potential generated by a homogeneous sphere has a com- 
pletely different behavior from that given in equation (4) , 
while ^p""" as well as $p^ compare better to a smoothed 
point-mass potential. This is clearly seen in Figure 1 and 
formally shown in § 2.1.2 and § 2.1.3. Thereupon, one 
should expect that computation results should depend 
only marginally on the adopted gravitational potential, 
except when #p = $ 



HS 



2.1.1. Homogeneous Sphere Solution 

The gravitational potential generated by a homoge- 
neous spherical distribution of matter can be calculated 
in a straightforward way by a direct integration of the 
gravitational force. Thereby, one finds 



GMp 
2 So 



if 5 < 5p 



iiS>Sp 



(6) 



where Sp is the radius of the sphere, i.e., the planet's 
radius. No smoothing is needed in this case since the force 



^Thc radius of the Hill sphere is Rjj = a where a is 

the distanec of the star from the planet. Considering only the 
linear terms in q, this length is equal to the distance from the 
inner Lagrangian point (LI) to the secondary. Although we 
named S the distance from the planet, we preferred to keep a 
more familiar notation to indicate the Hill radius. 



converges linearly to zero as the distance S approaches 
zero. Thus, there is no risk of numerical instabilities. 

Strictly speaking, equation (6) is valid inside very ex- 
tended and nearly homogeneous envelopes without consid- 
erable cores. Then, one may think of the functions $™ 
and <l?p^ as rendering two opposite extreme situations. 
Though 3>p ® does not represent a very realistic scenario 
of planet formation, for the sake of comparison and com- 
pleteness, we will apply this potential to high- as well as 
low-mass bodies. 

2.1.2. Stevenson's Solution 

Stevenson (1982) proposed a simplified analytical 
model of protoplanets having envelopes with constant 
opacity and surrounding an accreting solid core. Ho de- 
veloped a radiative zero solution for hydrostatic and fully 
radiative spherical envelopes, which implies that both hy- 
drostatic and thermal equilibrium are assumed inside the 
planet's atmosphere. Under these hypotheses, the core 
can grow up to a critical mass whose value is that beyond 
which at least one of the two equilibriums ceases to exist 
(see discussion in Wuchterl 1991) and the structure can- 
not be strictly static amy longer. The critical core mass 
also sets an upper limit to the envelope and total mass of 
the planet. It can be proved that, at this critical point, 
Mc/Mp = 3/4, where the mass of the core Mc is influenced 
by neither the nebula density nor its temperature. 

The potential of the gravitational field established by a 
fully radiative envelope can be obtained from the density 
profile (see Stevenson 1982) by applying the Poisson equa- 
tion. Since the solid core size is by far below the resolution 
limit of these computations, the form of the Stevenson's 
potential can be cast in the form 

GMe QMt 



(7) 



iiS>Sp 



In equation (7) we have indicated with Sc the core 
radius. The quantity is equal to the planet's envelope 
mass Me = Mp — Mc divided hy \n(Sp/ Sc). The presence 
of the parameter 5, also in the solution valid outside the 
envelope, is necessary for continuity reasons at 5 = Sp. 
In these simulations we set 5 = 0.05 -Rh. 

If the core has a density Pc = 5.5 gcm~^ and accretes 
at the rate of 5 x 10~^ Meyr"^, assuming an envelope 
opacity equal to 1 cm^ g"'^ and a mean molecular weight 
of 2.2, the critical total mass is 36 A/®. Hence we will use 
equation (7) only for protoplanets whose total mass Mp 
is less than that value. Furthermore, we will suppose that 
the ratio of the total planetary mass to the core mass is 
the critical one. Therefore the core mass is always known 
once the mass Mp is assigned a value. Then, supplying 
Pc the radius Sc can be fixed. 
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Fig. 1. - Left panel. Gravitational acceleration {—d^p/dS) inside a 20 planet as generated by the three 
potential functions ^p®, and ^p""" (eqs. [4], [6], and [7], respectively). Accelerations are normalized to 

—GMp/Sp, where Sp is the envelope radius. The core mass (see § 2.1.2) is 15 while Sp = 0.52 i?H- The core 
is supposed to have a density of 5.5 gcm~^. The point-mass potential is smoothed over a length s = 0.1 i?H 
whereas in the Stevenson's potential 6 = 5 x 10~^i?H- Right panel. The same quantity is displayed inside a 
90 M© planet but this time involving the potential solution $p^ (eq. [9]). The adiabatic exponent is F = 1.43 
while the core mass is 60 . In all circumstances, the gravitational potential outside the envelope radius is of 
the type given in equation (4). 



The effects caused by equation (7) arc actually similar 
to those caused by equation (4), as soon in Figure 1. In a 
more formal way, inside of the sphere S = Sp, the normal- 
ized difference TZ^^ between the two gravitational fields 
can be quantified by the ratio of \d^l'^/dS - d^l^/dS\ 
that is 



to d^^'^'/dS^ 



ST 



1 



4 \a[Sp/Sc) 

2 



+ ln 



S 



(8) 



In the above relation, the equality e = 5 was imposed. 
TZ^^ is a decreasing function of S. Referring to the models 
addressed in the left panel of Figure 1, 7^^'^ ~ 5% at 
S/S-p = 1/3. This result is only slightly affected by the 
mass of the planet because Sp is a slowly varying function 
of Mp. 

2.1.3. WuchterVs Solution 

Along with the fully radiative envelope, other static so- 
lutions were found. Following the track of Stevenson's ar- 
guments, Wuchterl (1993) developed an analytical model 
for protoplanets with spherically symmetric and fully con- 
vective envelopes. In this case the hydrostatic structure 
is determined by the constant entropy requirement that is 
appropriate when convection is very efficient. Supposing 



that the adiabatic exponent F = dlnp/dlnp is constant 
throughout the envelope, integrating the envelope density 
one finds that a solution for the planet gravitational po- 
tential is 



(AX*)- (AX*) 



n 



iiS<Sl 



lis > si 



(9) 

where we set C = (SF - 4)/(F - 1) and H = (Sc/S^)^. 
Moreover, the envelope mass is written as Me — (1 — 
n) M| . The condition for stability of gas spheres (F > 
4/3) implies that C, is positive. This particular form of 
is obtained by choosing an envelope radius equal to 
Si = (F - l)i?ac- The length iJac = GMp/cl is called 
"planetary accretion radius" . Outside the sphere S = i?ac 
the thermal energy of the gas is higher than the gravita- 
tional energy binding it to the planet. 

As in Stevenson's solution, critical mass values exist for 

the envelope structure to bo static. However, unlike the 
fully radiative envelope case, now the critical core mass 
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depends on both the temperature and the density of am- 
bient material. Furthermore, the critical mass ratio is 
Mc/Mp = 2/3 (for details see Wuchterl 1993). Setting T 
to 1.43 and the mean molecular weight to 2.2, if nebula 
conditions are TNob ~ 100 K and pNob = 10"^" gcm~^, 
the critical total mass is Mp = 274 M®. 

Wuchterl's solution well applies to massive protoplan- 
ets since they are likely to bear convective envelopes. 

Concerning the differences between equation (9) and 
equation (4), the situation is alike to that met in § 2.1.2 
(see right panel of Fig. 1). In this circumstance, for S < 
Sp, the normalized difference can be written as 



i-n 
(- 



c 



C-1 



C- 1 



(10) 



As before, the equality e = 5 is assumed and 71^^ de- 
creases with S. Using the parameters adopted for the 
models illustrated in the right panel of Figure 1, TZ^^ ~ 
18% at S/Sp = 1/3. In the considered range of masses, 
this number is nearly constant. In fact, equation (10) can 
be approximated to 



(11) 



2.2. Physical Parameters 

We consider a protostellar disk orbiting a one solar- 
mass star. The simulated region extends for 2 it around 
the polar axis and, radially, from 2.08 to 13 AU. The 
aspect ratio is fixed to h — 0.05 throughout these compu- 
tations. As in Paper I, the disk is assumed to be symmet- 
ric with respect to its midplane. This allows us to reduce 
the latitude range to the northern hemisphere only, where 
varies between 80° and 90°. The co-latitude interval 
includes 3.5 disk scale-heights and therefore it assures a 
vertical density drop of more than six orders of magni- 
tude (see § 3.3). The mass enclosed within this domain 
is Md = 3.5 X 10~^ Mq, which implies, in our case, that 
a mass equal to 0.01 Mq is confined inside 26 AU. Disk 
material is supposed to have a constant kinematic viscos- 
ity 1/ = 10^® cm^s~^, corresponding to a = 4 x 10~^ at 
the planet's location. 

The orbital radius of the planet is Rp = 5.2 AU and its 
azimuthal position is fixed to (p = ipp = ir . We concentrate 
on a mass range stretching from 1.5 M© to one Jupiter- 
mass (Mi^), implying that the mass ratio q = Mp/Mi, € 
[4.5 X 10-*',10-^], if M* = 1 Mq. a detailed list of 
the examined planetary masses, along with the adopted 
potential form is given in Table 1. 

The choice of few of the above parameters represents 
one typical example during the early phase of planet for- 
mation. Additionally, these simulations offer the good 




Fig. 2. — Envelope radius S'p, as provided by P. Bo- 
denheimer (2001, private communication), compared 
to the Hill radius i?H and the accretion radius i?acc- 
In a locally isothermal disk with aspect ratio h, Ru 
is larger than i?ac when q < /i^/\/3. While the scale 
of the left vertical axis is referred to orbital radius 
Rp = 5.2 AU, the scale on the right one is compared 
to the radius of Jupiter Rq^ = 7.1 x 10^ km. 



advantage that the system of equations is cast in a non- 
dimensional form, thus all of the outcomes are "scale-free" 
with respect to M^^, Md, and Rp. 

According to studies of the early evolution of proto- 
planets, the Roche lobe is usually filled during the growth 
phase. In such calculations the envelope is allowed to ox- 
tend to either the Hill radius Rn or the accretion radius 
Rac (Bodenheimer & Pollack 1986; Wuchterl 1991; Tajima 
& Nakagawa 1997). Except for Wuchterl's solution, where 
we set Sp = Sp = (F — 1) i?ac, the estimates of planet radii 
used in the simulations were provided by P. Bodenheimer 
(2001, private communication). They originate from a 
combination of Rn and /Jac at an ambient temperature of 
T = 100 K (see Fig. 2). The values of Sp, employed in 
each model, are also reported in Table 1. 

It is worthwhile to stress that the planetary (or more 
properly, the envelope) radius Sp does not represent any 
real physical boundary but only the distance beyond 
which the planet's potential reduces to the one given in 
equation (4), i.e., to a point-mass potential. 

3. Numerical Issues 

The set of hydrodynamic equations that characterizes 
the temporal evolution of a disk-planet system is solved 
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Table 1 

Planetary Masses and Adopted Gravitational Potential. 
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Sp/Ru 


Ru/a 




Potential 


333 


1.00 


X 


10- 


-3 


0.87 


6.9 X 


10- 


-2 


PM, HS 


253 


7.56 


X 


10- 


-4 


2.30*^ 


6.3 X 


10- 


-2 


KW 


166 


5.00 


X 


10- 


-4 


0.78, 1.70^ 


5.5 X 


10- 


-2 


PM, KW, HS 


93 


2.83 


X 


10' 


-4 


1.20'^ 


4.5 X 


10- 


-2 


KW 


67 


2.00 


X 


10- 


-4 


0.70, 0.96'^ 


4.0 X 


10- 


-2 


PM, KW 


33 


1.00 


X 


10- 


-4 


0.60 


3.2 X 


10- 


-2 


PM, HS 


29 


8.80 


X 


10- 


-5 


0.58 


3.1 X 


10- 


-2 


ST 


20 


6.00 


X 


10- 


-5 


0.52 


2.7 X 


10- 


-2 


HS, ST 


15 


4.50 


X 


10- 


-5 


0.46 


2.5 X 


10- 


-2 


ST 


12.5 


3.75 


X 


10- 


-5 


0.44 


2.3 X 


10- 


-2 


ST 


10 


3.00 


X 


10- 


-5 


0.38 


2.1 X 


10- 


-2 


PM, HS, ST 


7 


2.10 


X 


10- 


-5 


0.34 


1.9 X 


10- 


-2 


ST 


5 


1.50 


X 


IQ- 


-5 


0.29 


1.7 X 


10- 


-2 


HS, ST 


3 


1.00 


X 


10- 


-5 


0.23 


1.5 X 


10- 


-2 


ST 


1.5 


4.50 


X 


10- 


-6 


0.16 


1.1 X 


10- 


-2 


ST 



Values are rounded to the nearest integer numbers. 

^Planetary radius used in the Wuchterl's solution: Sp = = 
(r- l)i?ac (see § 9). 

Note. — List of all the simulated planet masses: q = Mp/M^ is 
the non-dimensional quantity that enters the simulation. Note that 
Mp = 333 Me = 1.05 M^. The ratio Rn/a is equal to 
(sec § 2.1). Unless stated otherwise, envelope radii Sp are expressed 
through a combination of the Hill (Rh) and the accretion (-Rac) radius 
of the planet, as plotted in Figure 2 (courtesy of P. Bodenheimer) . 
When required, we set Mc = (2/3) Mp for a fully convective envelope 
and Mc = (3/4) Mp for a fully radiative one. The core radius is 
computed assuming a constant density pc = 5.5 gcm~^ in both cases. 
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numerically by means of a finite difference scheme pro- 
vided by an early FORTRAN-coded version of Nirvana 
(Ziegler & Yorke 1997; Ziegler 1998). The code has been 
modified and adapted to our purposes as described in Pa- 
per I, Paper II, and references therein. In order to inves- 
tigate thoroughly the flow dynamics in the neighborhood 
of the planet, a sufficient numerical resolution is required. 
We accomplish that by employing a nested-grid strategy. 
This can bo pictured as either a set of grids, each hosting 
an inner one, or a pyramid of levels: the main grid (level 
one) includes the whole computational domain, while in- 
ner grids (higher levels) enclose smaller volumes around 
the planet, with increasing resolution. Levels greater than 
one are also called subgrids since they are usually smaller 
in size. The linear resolution, in each direction, doubles 
when passing from a grid to the inner one. 

The basic principles upon which the nested-grid tech- 
nique relies and how it is applied to disk-planet simu- 
lations in two dimensions is explained in detail in Pa- 
per II and references therein. The extension to the three- 
dimensional geometry, though requiring some more com- 
plexity in the exchange of information from one grid level 
to the neighboring ones, is nearly straightforward. Coarse- 
fine grid interaction (which sots the boundary conditions 
necessary for the integration of subgrids) is accomplished 
via a direction-splitting procedure. Thus, with respect 
to the discussion in Paper II, the addition of a third 
dimension just requires the implementation of the third 
direction-splitting step in the algorithm. Fine-coarse grid 
interaction (which upgrades the solution on finer resolu- 
tion regions) involves, in 3D, volume- weighted averages in 
a spherical polar topology. These are given in the Ap- 
pendix. 

3.1. General Setup 

For the study of the variety of planetary masses in- 
dicated in Table 1, meeting both the requirements of 
high resolution and affordable computing times, we re- 
alized a series of six grid hierarchies, whose character- 
istics are given in Table 2. With mass ratios q larger 
than 2 x 10"* (67 Me) only grids GO and Gl are utifized 
whereas smaller bodies are investigated with the other 
grid hierarchies. Thereupon, the finest resolution we ob- 
tain in the whole set of simulations varies form 0.03 to 
0.06 Rh- In all of the models presented here, the planet is 
centered at the corner of a main grid cell, which property 
is retained on any higher hosted subgrid. As the planet 
radial distance is _Rp = a/(l -I- q), we adjust it by tun- 
ing the value of the star-planet distance a. Adjustments 
never exceed 0.7% over the nominal value of Rp given in 
§ 2.2. Every model is evolved at least till 200 orbits. The 
evolution of massive planets (Mp > 67 M®) is followed 
for 300 to 400 orbital periods because they take longer to 
settle on a quasi-stationary state. 

Gas accretion is estimated following the procedure 
sketched in Paper II. For better accuracy, mass is re- 
moved only from the finest grid level according to an ac- 



cretion sphere radius Kac and an evacuation parameter 
Kev The former defines the spherical volume which con- 
tributes to the accretion process whereas the latter can 
be regarded as a measure of the removal time scale in 
such volume. Two-dimensional simulations showed that 
the procedure is fairly stable against these two parame- 
ters. We constrain the amount of removed mass per unit 
volume not to exceed 1% of that available, as was done in 
Paper II. Regarding the extension of the sphere of accre- 
tion /vac, we performed simulations using different values, 
as stated in Table 3. Since the planet actually works as a 
sink, our procedure only furnishes upper limits to realis- 
tic planetary accretion rates (see discussion in Tanigawa 
& Watanabe 2002). 

However, we also inquire how mass removal can pos- 
sibly affect gravitational torques and, more generally, the 
dynamics of the flow in the planet neighborhood by means 
of models in which accretion is prohibited. 

3.2. Boundary Conditions 

In order to mimic the accretion of the flow towards the 
central star, an outflow boundary condition is applied at 
the inner radial border of the computational domain. The 
outer radial border is closed, i.e., no material can flow in 
or out of it. The same condition exists at the highest lati- 
tude 9 = 80° . Since the disk is symmetric with respect to 
its midplane as mentioned in § 2.2, symmetry conditions 
are set at 5 = 90°. On subgrids, except for the midplane 
where symmetry conditions are applied, boundary values 
are interpolated from hosting grids, by means of a mono- 
tonised second-order algorithm (see Paper II for details). 

The open inner radial boundary causes the disk to 
slowly deplete during its evolution. For all the models 
under study, we observe a depletion rate Md = — M* 
measuring w 10~^ M0yr~^, in agreement with the ex- 
pectations of stationary Keplerian disks: M*^ = Sttz/E 
(Lynden-Bell & Pringle 1974). 

In cases of gap formation, material residing inside the 

planet's orbit tends to drain out of the computational do- 
main. Since this material transfers angular momentum to 
the planet, the lack of it may contribute to reduce both 
the migration time scale and the planet's accretion rate. 
To evaluate these effects, a Jupiter-mass model was run 
with a closed (i.e., reflective) inner radial border. 

3.3. Initial Conditions 

The initial density distribution is a power-law of the 
distance from the rotational axis r = R sin times a Gaus- 
sian in the vertical direction 



p{t = 0) = po{r) exp 



cot ( 



(12) 



which is appropriate for a thin disk in thermal and hydro- 
static vertical equilibrium. The dependency of the mid- 
plane value po with respect to r is such that the initial 
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Table 2 

Grid Hierarchies Utilized in the Simulations. 



Grid 


Main Grid Size 


ng 


Subgrid Size 


No. of models 


GO 


121 X 13 X 319 


4 


54 X 12 X 48 


10 


Gl 


143 X 13 X 423 


4 


64 X 12 X 64 


5 


G2 


121 X 13 X 319 


5 


54 X 12 X 48 


19 


G3 


143 X 13 X 423 


5 


64 X 12 X 64 


12 


G4 


121 X 23 X 319 


5 


54 X 22 X 48 


4 


G5 


133 X 13 X 395 


5 


84 X 16 X 84 


1 



Note. Grid sizes are reported as the number of grid points per 
direction: xNgX N^. The third cohimn {ng) indicates the num- 
ber of levels within the hierarchy. In order to achieve sufficient res- 
olution within the Roche lobe of the planet, grids GO and Gl have 
been used only for planetary masses in the range [67 , 1 Mq^] . 
Grids are ordered according to their computing time requirements, 
which grow from top to bottom. The hierarchy G5 has been em- 
ployed to execute the model with Mp = 12.5 (see discussion at 
the end of § 4.2). 



surface density profile E decays as l/\/r, as required by 
the constant liinematic viscosity. However, we also ran 
a few models where the relation pair) is such to account 
for an axisymmctric gap, as often done to speed up the 
computations at early evofutionary times. 

The initial velocity field of the fluid is a purely, counter- 
clockwise, Keplerian one corrected by the grid rotation: 
u{t = 0) = {0,0, VK — Opr). Thus, the partial support 
due to the radial pressure gradient is neglected in the be- 
ginning. 

4. Simulation Results 

4.1. Flow Dynamics near Protoplanets 

Two-dimensional computations have shown that a cir- 
cumplanetary disk forms around Jupiter-type planets, 
extending over the size of its Roche lobe (Kley 1999; 
Lubow, Seibert, & Artymowicz 1999; Tanigawa & Watan- 
abe 2002). The numerical experiments conducted in PA- 
PER II proved this characteristic to belong not only to 
massive bodies but also to protoplanets as small as 3 . 
The authors demonstrated that the flow of such disks is 
approximately Keplerian down to distances 0.1 -Rh from 
the planet. One of the main features of these disks is a 
two-arm spiral shock wave whose opening angle (that be- 
tween the wave front and the direction toward the planet) 
is an increasing function of the mass ratio q and, below 
Mp = 67 A/0, it is roughly given by arctan(7V/f), in which 
A4 = |m|/cs is the Mach number of the circumplanetary 



flow (D'Angelo, Kley, & Henning 2002). The spiral pat- 
terns shorten and straighten as the perturber mass de- 
creases. Eventually, for even smaller masses they disap- 
pear and are not observable anymore when one Earth- 
mass is reached. 

Three-dimensional simulations shed new light on 
these circumplanetary disks, demonstrating that they 
can behave somehow differently from what depicted by 
two-dimensional descriptions. Differences become more 
marked as the mass of the embedded protoplanet is low- 
ered. A major point is that spiral waves are not so pre- 
dominant as they arc in the 2D geometry. This is clearly 
seen in the top rows of Figure 3 and 4, where the midplane 
{$ = 7r/2) density is displayed for four different planetary 
masses. The double pattern of the spiral is still visible 
around a 67 protoplanet (Fig. 3, top right panel). 
However, when considering a planet half of that size, spi- 
ral traces are too feeble to be seen on the image (Fig. 4, 
top left panel). Such an occurrence was to be expected 
since the energy of the flow is not only converted into the 
equatorial motion but can be also transferred to the ver- 
tical motion of the fluid. It was already known from wave 
theories for circumstellar disks (Lubow 1981; Lubow & 
Pringle 1993; Ogilvie & Lubow 1999) and related numeri- 
cal calculations (Makita, Miyawaki, & Matsuda 2000) that 
the three-dimensional propagation of spiral perturbations 
may be significantly different from that obtained in two 
dimensions because of the existence of vertical resonances. 
Furthermore, Miyoshi et al. (1999) already noticed in their 
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Table 3 
Mass Accretion Parameter k-, 







q 






Accreting Only^ 


333 


1.00 


X 


10-3 


0.20, 0.15, 0.10 


No 


253 


7.60 


X 


10-4 


0.10 


Yes 


166 


5.00 


X 


10-4 


0.10 


Yes 


93 


2.80 


X 


10-4 


0.10 


Yes 


67 


2.00 


X 


10-4 


0.20, 0.10 


No 


33 


1.00 


X 


10-4 


0.20, 0.10 


No 


29 


8.80 


X 


10-5 


0.20, 0.10 


No 


20 


6.00 


X 


10-5 


0.15, 0.10 


No 


15 


4.50 


X 


10-5 


0.10 


Yes 


12.5 


3.75 


X 


10-5 


0.10 


Yes 


10 


3.00 


X 


10-5 


0.10 


No 


7 


2.10 


X 


10-5 


0.10 


Yes 


5 


1.50 


X 


10-5 


0.10 


No 


3 


1.00 


X 


10-5 


0.10 


Yes 


1.5 


4.50 


X 


10-*^ 


0.07 


Yes 



'^"No" entry stands for the existence of a non-accreting 
model. 

Note. — The parameter represents the radius of the ac- 
creting region. Within this sphere the mass density is reduced 
by roughly 1% after every time step. The length k^c = 0.1 i?H 
should be small enough to make the accretion procedure al- 
most independent of the evacuation parameter Kov (Tanigawa 
k Watanabe 2002). At « = 3 x 10-5 (Mp = 10 M©) a non- 
accreting simulation was performed with $p = (f>p^ (eq. [6]) as 
well as with $p = $p'^ (eq. [7]). In the case of lowest mass 
model (Mp = 1.5 M®), we allowed the ratio Kac/Sp to be less 
than 0.5, as in all of the other models. For a better evaluation 
of Mp, we used a modified version of the grid system G3 which 
contains a sixth level, comprising (approximately) the planet's 
Hill sphere. 
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Fig. 3. — Density slices with an overplotted two-coinponcnt velocity field. In order to comprise the whole 
density range, the color scale represents the logarithm of the density (logp). The evolutionary time is 200 
orbital periods. Slices are cut at the planet location: 9 = 7r/2 (top), (p — (pp (middle), and R = Rp (bottom). 
Both models are accreting and executed with Wuchtcrl's potential. In physical units, p = 10~^ is equivalent to 
4.22 X 10~^^ gcm""^. Left panels. Planet with a total mass Mp = 166 and a core mass equal to 2/3 Mp 
(critical core mass). Right panels. Planet with Mp — 67 and a critical core mass. The curve drawn in 
each panel indicates the trace of the Roche lobe. In the left and right top panels, the flow reaches velocities equal 
to ~ 3 and ~ 2 kms~^, respectively. In the other panels, maximum velocities are on the order of ~ 4 kms~^. 
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shearing sheet models the weakening effect of the finite 
disk thickness upon the formation of spiral waves around 
embedded protoplanets. Hence, we may argue that the 
averaging of the pressure and the gravitational poten- 
tial, which is accomplished in an infinitesimally thin disk, 
enhances spiral features in disks. 

The remainder of this section is devoted to a general 
description of the vertical circulation of the material in 
the vicinity of the planet. We start inspecting what hap- 
pens in the slice ip = ipp, i.e., in the R-9 surface contain- 
ing the planet (see middle panels in Fig. 3 and 4). The 
first thing to note is that the material above the equa- 
torial plane moves toward the planet with a negligible 
meridional component, in fact -C u (where u = 
Therefore matter is nearly confined to ^-constant surfaces. 
This circumstance favors the use of 2D outcomes as pre- 
dictions of 3D expectations (Masset 2002). However this 
turns out to be true only far away from the planet. In 
fact, as the fluid enters a certain region around the Hill 
sphere, its dynamics changes drastically. The beginning 
of this zone is marked by two shock fronts, which actually 
develop well outside the Roche lobe of the restricted three 
body problem. The distance of the shock fronts from the 
perturber, if compared to Rh, shrinks from 2.61 to 1.16 as 
Mp is increased from 10 to 1 Mi^. Generally, shocks 
are not placed symmetrically with respect to the planet. 
Past these shocks, material is doficctcd upward and then 
it recirculates downward, while approaching the radial po- 
sition of the planet (see middle panels in Fig. 3 and 4). 

At 7? « -Rp, matter suffers an unbalanced gravitational 
attraction by the planet and accelerates downward to- 
wards it. Velocities are supersonic, reaching a Mach num- 
ber ~ 8, when 67 M® < Mp < 1 Mi^, and ~ 2 
if Mp = 20 M®. They become subsonic for planetary 
masses between 10 and 1.5 M®, ranging from 30% to 5% 
of the local sound speed. Because of the sinking material, 
at S = 7r/2 the fiow field is slightly horizontally divergent 
from the planet location. 

As one can judge from Figure 4 (middle panels), the 
flow gets less and less symmetric, with respect to Rp, as 
the planet mass is reduced. Recirculation persists before 
the planet {R < Rp) but vanishes behind it {R > Rp). 
Any symmetry disappears starting from Mp = 5 M® 
downwards. 

Along the azimuthal direction (slice R = Rp, i.e., in 

the surface Lp 9 containing the planet), the bottom pan- 
els of Figure 3 and 4 show an even more complex situa- 
tion. Below 166 M®, the region of influence of the planet 
appeaxs to be more comparable in size with the Roche 
lobe. However, apart from that, the general behavior of 
the flow differs from case to case, having iu common a 
rapid descending motion when the material lies above the 
planet. Around Jupiter-sized planets some kind of weak, 
non-closed, recirculation may be seen. This flow feature 
is still present at both sides of a 29 M® planet (bottom 
left panel. Fig. 4), whereas it tends to vanish in models 
with lower g-ratios. 



4.1.1. N on- accreting Protoplanets 

Here we should dedicate some attention to the differ- 
ences existing between accreting and non-accrotiug \iio- 
toplanets. Since the gas is locally isothermal, pressure 
is proportional to the density according to equation (2). 
Because of the mass removal, density nearby the planet 
is lower in accreting models than it is in non-accreting 
ones. In Table 4 the mass A#e enclosed within the enve- 
lope radius Sp is quoted for the two sets of models, along 
with the mean density po. These values demonstrate that 
the amount of material contained in the volumes of non- 
accreting planets can be considerably larger than in the 
other case (even 6 times as much). Since the pressure 
must converge in the two cases, when the distance from 
the planet 5 S> Sp, a higher mean pressure in the enve- 
lope intuitively implies a larger magnitude of the pressure 
gradient inside this region. 

As an example, we illustrate in Figure 5 the velocity 
field in two perpendicular slices 6 = 1^/2 (i.e., the equar 
torial plane) and ip = ipp (i.e., the surface R-6 passing 
through the planet), in order to show how the enhanced 
density values affect the local circulation. The targeted 
body has Mp = 20 M® because, among the eight avail- 
able models for which accretion is not considered (see Ta- 
ble 3), this is the one that suffers an orbital migration 
substantially different from accreting counterpart mod- 
els. From the isodensity lines displayed in Figure 5, one 
can infer that matter is spherically distributed around the 
non-accreting protoplanet (left panels). For this reason 
the net torque arising within a region of radius w Rn is 
nearly zero. This docs not happen in the other case be- 
cause the symmetry is not so strict. 

The upper right panel clearly indicates the existence 
of a rough balance between gravitational and centrifugal 
force, with the pressure gradient playing a marginal role 
in opposing the planet potential gradient. On the other 
hand, from the circulation in the upper left panel of Fig- 
ure 5 one can deduce that the pressure gradient is no 
longer negligible compared to the potential gradient and 
can therefore counterbalance its effects. 

Moreover, the flow above the disk midplane (center and 
lower left panels) suggests that gas is ejected at R < Rp. 
Such phenomenon must be ascribed to the pressure gra- 
dient as well, since the fluid opposes any further compres- 
sion. These qualitative arguments will be quantitatively 
corroborated in § 5, where we will show that the increased 
amount of matter causes the envelope to be pressure sup- 
ported. 

We mention in the caption of Figure 5 that the flow 

may travel supersonically within the atmospheric region, 
though both Stevenson's and Wuchterl's gravitational po- 
tential rely on the hypothesis of quasi-hydrostatic equilib- 
rium. Such discrepancy can be attributed to the rate of 
mass removal from the innermost parts of the planet's en- 
velope, which is not considered in the derivation of those 
analytic solutions. In fact, while supersonic speeds have 
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Fig. 4. — Density slices of the same type as in Figure 3, illustrating low-mass accreting planets after 200 orbits. 
In these simulations, Stevenson's potential is employed. In physical units, the density scale is as in Figure 3. 
The density value p = 10~^ corresponds to 4.22 x 10~^^ gcm~^. Left panels. Planet with Mp = 29 and 
a critical core mass Mc ~ 21.7 M^. Maximum flow speeds, within S = i?H, are ~ 1 kms"""^ in each of the three 
panels. Right pcinels. Planet with Mp = 10 and a core mass and Mc = 7.5 Mq. Inside of the Hill sphere, 
maximum velocities are on the order of 0.2 kms~^. 
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Fig. 5. Two-component velocity field on the equatorial plane (top row ) and on the slice Lp = ip-p (middle 
and bottom rows) for a non-accreting, 20 planet (left panels), and an accreting planet of the same mass 
(right panels) with an accreting sphere radius Kac = 0.15 i?H- Lines of equal logp are also drawn. The dashed 
line represents the Roche lobe of the restricted three-body problem. In the accreting model, the equatorial flow 
{0 = 7r/2) in the planet's envelope has a Mach number M. < 1.2 whereas, for the meridional flow, < 2. 
In the non-accreting computation the fluid travels subsonically in both the midplane and the meridional slice 
{M < 0.2). 
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Table 4 

Mass enclosed within the Envelope Radius. 



No Accretion Accretion 





Me 














67 


7.07 X 10-1 


9.37 X 10- 


ll 


1.07 X 10" 


-1 




7.67 X 10-2 


29 


2.74 X 10-1 


1.45 X 10" 


■10 


6.49 X 10" 


-2 






20 


2.26 X 10-1 


2.43 X 10" 


-10 


3.68 X 10- 


-2 


3.18 X 10-2 




10 


1.29 X 10-=^ 


7.12 X 10" 


-11 


9.76 X 10- 


-3 






5 


2.24 X 10-3 


5.56 X lo- 


-11 


2.17 X 10- 


-3 







'^From models with = 0.10 i?H- 
'^From models with = 0.15 i?H- 
•^From models with = 0.20 i?H- 



Note. — All masses are relative to the Earth's mass. The mean density pe within 
the planet's radius is expressed in cgs units. We consider only simulations which were 
run with the Stevenson's potential. The mass Mg is evaluated assuming a disk mass 
Md — 3-5 X 10^'^ Mq. Hence, the values in the Table do not account for the disk 
depletion rate Mdi which is on the order of 3 x 10-^ yr-^. 



been measured in some accreting models, the flow is al- 
ways subsonic in the envelope of non-accreting planets. 
This can be also understood with simple arguments. We 
said before in this section that mass accretion is responsi- 
ble for the Keplerian-like circulation around protoplanets, 
in the disk midplane. Hence, the local (equatorial) Mach 
number should be on the order of y/qR^JS/h. If evalu- 
ated a.t S ~ Sp, this quantity yields « 1.3 for a 20 
body, which is comparable to the value reported in the 
caption of Figure 5. Along the vertical direction, if the 
velocity is approximated to that of a spherically accreting 
flow then ue ~ Mp/iAivS^ p). At 5 = Kac = 0.1 Rn, that 
relation gives a (meridional) Mach number M ~ 1.8 when 
applied to the accreting model shown in Figure 5. Also 
this number is similar to the measured value. 

4.2. Gravitational Torques and Orbital Mi- 
gration 

Gravitational torques are believed to be responsible for 
the migration of protoplanets from their initial formation 
sites. In this work torques arc directly estimated from 
the gravitational force exerted by each fluid element of 
the circumstellar disk on the planet. When computing 
the gravitational force, we consider the density solution 
on the finest available grid level. This procedure permits 
to obtain higher accuracy because of the increasing reso- 
lution of hierarchy levels. 



Since the planet is an extended object, torques act- 
ing on each of its portions should be calculated and then 
added up to give the total torque vector, whose most gen- 
eral expression® is 



(Jip + S) 

GdMujR) dMpjS) 
\R-Rp-S\^ 



R. 



(13) 



Equation (13) explicitly states that circumstellar mate- 
rial can alter both the orbital and rotational spins of a 
protoplanet. However, here we shall confine our study 
to variations of the planet's orbital angular momentum 
because the evaluation of the rotational spin requires a 
rigorous treatment of the envelope self-gravity. This is 
not done here, as stated in § 2 (see eq. [1]). Thereby, we 
can proceed as if the whole planetary mass were concen- 
trated in its geometrical center R = Rp and integrate the 
force over the whole disk domain, excluding the volume 
Vp occupied by the planetary envelope. Yet, outside of 
such volume the gravitational potential is always that of 
a point-mass object (see the behavior of eqs. [4], [6], [7], 



^Wc would like to stress that in this form all of the gravita- 
tional contributions (due to star, the planet, disk self-gravity, 
etc.) are implicitly enclosed in the differentials dMD{R) and 
dMp(S). 
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and [9], for S > Sp), therefore equation (13) simplifies and 

becomes 

Td = -RpX / V$™p(i?)dy(i?). (14) 

The orbital angular momentum of a protoplanet can 
be affected only by the 2:-component (that parallel to the 
polar axis) of the torque vector To. To avoid useless dis- 
tinctions, wo indicate this component as Tu. The sign 
of To determines the gain (positive torques) or loss (neg- 
ative torques) of orbital spin. Larger spins correspond 
to more distant orbits. Since torques generally change on 
time scales on the order of ~ 50 revolutions, we work with 
their final magnitudes. We usually observe a slow decay 
of |7d| with time (see also the end of § 4.3). 

Two-dimensional computations reveal that torques ex- 
erted by circumplanetary material may amount to a fair 
fraction of the total torque, unless a suitable smoothing 
length (usually of the size of the Hill radius) is used in the 
planet gravitational potential $p. This is caused by the 
high surface densities reached around the planet and the 
lack of a vertical torque decay which naturally occurs in 
three dimensions (see § 5). In fact, in 3D, we observe that 
torques arising from locations close to the planet do not 
play such an important role as they do in 2D. 

Analyzing the relative strength of torques exercised by 
different disk portions, it turns out that in the mass range 
q € [2x W~*, 1 X 10~^], dominating negative torques arise 
from distances S > 1.2 hRp, where h is the disk aspect 
ratio. Below 33 M®, the largest net contributions are 
generated by material lying between S ~ 0.6 h Rp and 
S ~ 2.2 hRp. Therefore we can conclude that predomi- 
nant torques are exerted at distances from the planet com- 
parable with the Hill radius. Not more than 10% of Tu is 
built up by matter located within S w 0.6 hRp. 

Apart from the 20 protoplanet, the total torque 
evaluated in non-accreting models does not deviate con- 
siderably from that estimated in accreting ones, indepen- 
dently of the used potential. In fact, simulations based 
on the potentials $p^, and (eqs. [4], [6], and 

[7], respectively) supply values of To which differ by less 
than 40%. This circumstance may signify that, whether 
or not a protoplanet is still accreting matter from its 
surroundings, this is not generally crucial to the gravi- 
tational torques by the circumstellar disk. Thereby, be- 
ing an exception, the case Mp = 20 deserves some 
comments. For such mass, the torque integrated over 
the first two grid levels (S > 2.2 hRp) yields a positive 
value for both accreting and non-accreting planets. When 
adding the contributions from the third and forth level 
(0.6/i-Rp < 5 < 2.2 HRp), the torque experienced by the 
planet lowers but, while it becomes negative in the ac- 
creting case, it still remains positive in the non-accreting 
counterpart. It is especially matter residing between 1 
and 2Rh from the planet that builds up the difference. 
As material in the uppermost grid level {S < 0.6 hRp) 
of the non-accreting simulation does not exert any signif- 



icant torques (some little negative contribution is indeed 
measured in the accreting model). To keeps the positive 
sign although, in magnitude, it is eleven times as small 
as that evaluated in the accreting case. This phenomenon 
of torque reversal for non-accreting planets with masses 
of about 20 Me may be related to the very long migra- 
tion time scales obtained for fully accreting models having 
masses Mp w 10 (see below, and Fig. 6). 

Conservation of orbital angular momentum implies 
that a protoplanet has to adjust its orbital distance from 
the central star because of external torques exerted by 
the disk. If the orbit remains circular, the time scale over 
which this radial drift motion happens is inversely pro- 
portional to To , according to the formula: 

a Mp a'^ Qp 

In equation (15) we indicated with f2p and a the planet's 
angular velocity and its distance from the star, respec- 
tively. Linear, analytical theories (e.g., Ward 1997) pro- 
vide two separate regimes governing the migration of low- 
{type I) and high-mass {type II) protoplanets. Both mi- 
gration types predict that the planet moves toward the 
star. More recent studies by Masset (2001) and Tanaka, 
Takeuchi, & Ward (2002) have reconsidered the role of 
co-orbital corotation torques and proved that they can 
be very effective in opposing Lindblad torques. Hence, 
they can significantly slow down inward migration. Two- 
dimensional results presented in Paper II well fit to these 
predictions. A further reduction of the migration speed 
is expected from a full 3D treatment of torques, as also 
derived numerically by Miyoshi et al. (1999) and theoret- 
ically predicted by Tanaka, Takeuchi, & Ward (2002). 

Figure 6 shows our estimates for the migration time 
scale tm as computed for models having different masses 
and in which the potential solutions 3>™, 3'p^, and ^p"^ 
were adopted. We compare these values with the two 
analytical theories developed by Ward (1997) (solid line) 
and Tanaka, Takeuchi, & Ward (2002) (dashed fine). The 
first of them comprises both migration regimes, though ac- 
counting only for Lindblad torques. The second theory is 
fimited to type I migration, albeit it treats both Lindblad 
and corotation torques. Moreover, the first is explicitly 
two-dimensional whereas the second is applicable in two 
as well as three dimensions. 

As one can see from Figure 6, numerical results are very 
similar for Mp > 67 M^, yielding tm ~ 5 x 10^ years, 
whatever of the four gravitational potential is used (see 
also Table 5). While this time scale is consistent with 
Ward's 1997 description if Mp = 67 M®, for more mas- 
sive planets it is nearly two times as short as the theo- 
retical prediction. The depletion of the disk inside the 
planet's orbit is probably responsible for part of the dis- 
crepancy (see § 4.5), because Ward's theory assumes a 
disk with a constant unperturbed surface density. In the 
type I regime, our numerical experiments with $p — $p'^ 
(eq. [7]) well reproduce the behavior of the analytical 
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Fig. 6. — Migration time scale versus the mass ratio q = M-p/M^. Outcomes from simulations carried out with 
the various forms of gravitational potential are marked with different symbols. A.t q = 10~^ and 6 x 10~^, 
calculations carried out with the point-mass potential were initiated with a density gap (sec § 3.3). To avoid 
confusion, migration rates obtained from models with $p = (eq. [6]) are quoted in Table 5. The solid line 
represents the theoretical prediction by Ward (1997), which does not include corotation torques and 3D effects. 
Both are indeed accounted for in the analytical model by Tanaka, Takeuchi, & Ward (2002) (dashed line). The 
parameter it a? Sp/MD, needed to draw the analytic curves, was retrieved from the initial surface density profile 
(see § 3.3). The scale on the right vertical axis gives the migration rates of a protoplanet orbiting at 1 AU from 
the primary. 



curve when Mp = 30, 20, 5, 3.3, and 1.5 Me. 

Computations executed witii $p — (eq. [6]) prob- 
ably underestimate the magnitude of differential torques 
because of the much weaker gravitational field. Neverthe- 
less, for Jupiter-mass and Earth-mass protoplanets, mi- 
gration times yielded by these models well compare to 
those displayed in Figure 6, as proved by the values re- 
ported in Table 5. 

Significant deviations from the linear estimate of 
Tanaka, Takeuchi, & Ward (2002) are observed in the 
mass interval [7, 15] M®, where the migration time is 
longest at 10 M®. For this planet tm, estimated with 
Stevenson's as well as the point-mass potential, is thirty 



times as long as the theoretical description by Tanaka, 
Takeuchi, & Ward (2002) predicts. This depends on the 
strong positive torques arising at 5 > 2 ft _Rp which are 
not efficiently contrasted by negative ones generated in- 
side S ~ hRp. However, for this particular planetary 
mass, we obtain discrepant estimates from computations 
performed with different resolutions. In fact, the simula- 
tion carried out with the grid G2 yields a positive total 
torque acting on the planet, i.e., it predicts an outward 
migration, whereas models based on the higher resolution 
hierarchies G3 and G4 provide a negative total torque. 
Yet, the absolute value of Tb evaluated with grid G4 is a 
seventh of that achieved with grid G3. 
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Table 5 

Migration Times obtained with the Homogeneous Sphere Potential. 





q 






[years] 


Accretion 


No Accretion 


333 


1.0 X 10" 


-3 


5.04 X 10" 


4.80 X 10^ 


166 


5.0 X 10- 


-4 


4.76 X 10^ 




33 


1.0 X 10" 


-4 


2.10 X 10^ 


2.54 X 10^ 


20 


6.0 X 10" 


-5 


5.74 X lO-"^ 




10 


3.0 X 10- 


-5 


5.26 X 10^ 


4.78 X 10^ 


5 


1.5 X 10" 


-5 


6.19 X 10^ 





Note. Migration rates from accreting and non- 
accreting models in which $p = (eq. [6]). Time scales 
are sensibly different from those indicated in Figure 6 only 
at Mp = 33, 20, and 10 M®. 



Since we believe that gravitational torques are ac- 
counted for in a more accurate fashion by hierarchy G3 
than by grid G4 because of the arguments in § 4.5, we 
may rely more on the outcome of the hierarchy G3 (shown 
Fig. 6) rather than on the other two. We note that such 
migration time {tm = 3.3 x 10® years) is roughly the dou- 
ble of that supplied by the Mp = 10 M®, non-accreting 
model. 

We have further inquired into the matter by running a 
simulation with Mp = 12.5 Me and $p = $p^ (see Ta- 
ble 1). Based on the experience acquired with models exe- 
cuted with grids G3 and G4, we set up the high-resolution 
hierarchy G5 (see Table 2). This grid is designed to better 
resolve those regions responsible for the strongest torque 
contributions in the mass interval [7, 15] M©. As seen in 
Figure 6, the resulting migration rate follows the trend es- 
tablished by the other assessments in this range of masses. 

The interesting property that computed migration time 
scales are very long for 10 M® planets may be caused by 
non-linearity effects. We note, in fact, that in numerical 
simulations the first traces of a trough in the density struc- 
ture is observed around the same value of Mp and gap for- 
mation starts when disk-planet interactions become non- 
linear. However, this issue needs to be addressed more 
thoroughly with future computations. 

4.3. Mass Accretion 

Three-dimensional computations of one Jupiter-mass 
bodies provide estimates of the mass accretion rate Mp 
on the same order of magnitude as those obtained by two- 
dimensional ones (see Paper I). Two-dimensional calcula- 
tions performed by the authors reveal a mgiximum of the 



accretion rate, as function of the mass, around 0.5 M-}^ 
(see Paper 11). Yet, those estimates appear surprisingly 
high in the very low mass limit. Part of the reasons may 
lie in the assumed flat geometry which cannot account for 
the vertical density stratification. The present simulations 
overcome this restriction, hence they permit to evaluate 
also the effects due to the disk thickness. 

The values of Mp is plotted against the planetary mass 
in Figure 7. As comparison, estimates relative to models 
with different gravitational potential solutions are shown. 
The overall behavior of the data points resembles that re- 
ported in Paper II, with a peak around 0.3 Af^j.. For 
Mp = 1 Mi^_ the agreement between two and three- 
dimensional models is very good and not much discrep- 
ancy is seen down to Mp = 20 M® , since values are com- 
parable within a factor 3 (see § 4.4). Below this mass, 
however, the accretion rate rapidly declines, which drop is 
not observed in 2D outcomes. In fact, one can infer from 
Figure 7 that the dynamical range of Mp stretches for 
more than two orders of magnitude. By using model re- 
sults obtained applying the point-mass, Stevenson's, and 
Wuchterl's potential (eqs. [4], [7], and [9], respectively) 
the following approximate relation can be found: 



log 



M„ 



bo + bi log q +b2 (log q) 



(16) 



whose cocflicients are bo = -18.47 ± 0.76, 6i = -9.25 ± 
0.38, and 62 = -1.266 ± 0.046. Equation (16) holds as 
long as the mass ratio q € [4.5 x 10~®, 10~^] or, for a one 
solar-mass star, when 1.5 M® < Mp < 1 Mq^. Such an 
equation can be applied to scenarios studying the global 
long-term evolution of young planets. 
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Fig. 7. — Planet's accretion rate as function of the normalized planet mass q. Different symbols stand for the 
different forms of gravitational potential $p adopted in the computations. Apart from those models run with 
the homogeneous sphere potential (see eq. [6]), mass is removed from a volume, centered on the planet, 
with radius Kac = 0.1 Ru (see Table 3 for some details concerning the simulation with Mp = 1.5 Mq). Models 
for which $p = have Kac equal to 0.2 Ra if Mp > 20 M© and to 0.15 i?H if Mp = 20 M©; otherwise /tac is 
set to 0.1 i?H- 



Calculations in which the homogeneous sphere poten- 
tial (eq. [6]) is adopted yield accretion rates substan- 
tially lower (from 3 to 15 times) than those achieved when 
the other potential forms are employed. This is due to the 
weak gravitational attraction this potential exerts within 
the planet's envelope. As proved by Figure 1, the gravita- 
tional field can be 100 times as small as that established 
by the other three potential functions for S < Sp. Also 
in this circumstance, a relation similar to equation (16) 
exists for which the coefficients are bo = —19.42 ± 2.68, 
bi = -9.96 ± 1.41, and 62 = -1.42 ± 0.18. 

While the accretion rate is fairly stable with time for 
masses below 30 M®, it keeps reducing for higher masses. 
Between 67 and 0.8 Mq^, Mp drops by 10 to 20% 
during the last 50 orbits of the simulations. This is an 
indication of a deepening gap and a depleting disk. As 



for the dependency upon the accretion volume, from our 
numerical experiments it is found that doubling the radius 

Kac, the accretion rate grows at most by 30%. The smaller 
the planet mass, the less sensitive Mp is to the parameter 

4.4. Comparison with 2D Models 

In this Section we aim at comparing the migration time 
scale as well as the planet's accretion rate obtained in this 
work with those presented in Paper II. However, while 
two-dimensional estimates of Mp (Fig. 25 in PAPER II) 
are directly comparable to those plotted hero in Figure 7, 
the time scales tm shown in Figure 20 of Paper II are not 
completely consistent with those in Figure 6. Therefore, 
they need to be corrected. 

This is because in the present study torques are in- 
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Fig. 8. — Left panel. Migration rates as evaluated in two-dimensional (Paper II) and three-dimensional 
models (this paper). Right panel. Using the same sources as in the left panel, a comparison of the growth 
time scales tq = Mp/Mp, between 2D and 3D outcomes, is shown. Filled triangles, in both panels, indicate 
results obtained from three-dimensional models in which the planetary gravitational potential is $p^, or 
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tegrated all over the disk domain excluding the planet 
envelope, i.e., the sphere of radius S = Sp. Instead, in 

Paper II the excluded region has a radius « 0.1 -Rh (for 
details see Paper II, § 5.4). From Table 1 one can see 
that, above 33 M^, the envelope radius Sp can be much 
larger than a tenth of the Hill radius. 

Figure 8 illustrates the migration rate tm (left panel) 
and the growth time scale tq = Mp/Mp (right panel) as 
computed in the two geometries. Orbital migration esti- 
mated by means of 3D simulations is slower than that eval- 
uated in 2D calculations only below 33 {q = 10"'*). 
As for planet's accretion, the most important difference 
is the rapid drop, for Mp < 10 M®, observed in disks 
with thickness. The larger values of three-dimensional es- 
timates, measured in the range 10 < Mp < 1 Mq^_, 
are due to the gap which is not so deep as it is in two- 
dimensional models, hence the average density around the 
planet is higher. This occurrence can be partly attributed 
to gravitational potential effects that, as we mentioned in 
§ 4.1, are intensified by the fiat geometry approximation. 

4.5. Numerical Effects 

In Paper II it was found that, upon increasing the 
smoothing parameter, there is a reduction of the torques' 
mismatch, over a region around the planet whose lin- 
ear size is comparable with the double of the smoothing 
length. A 33 model was run with a point-mass po- 
tential without any kind of softening. This is possible be- 
cause none of the hydrodynamical variables is placed at a 



cell corner, where the planet dwells. A similar simulation 
was performed applying a grid dependent smoothing of 

the type described in Paper II. Resulting migration time 
scale and mass accretion are not significantly affected by 
the smoothing choice. 

As for the consequences of the circumstellar disk de- 
pletion, inside of the planet's orbit, we ran a Jupiter-mass 
model in which both inner and outer radial borders were 
closed. Since more material is available in the disk portion 
R < Rp (roughly twelve times as much), one should ex- 
pect larger values for both Mp and tm- Indeed, accretion 
is two times as much as that calculated in the model with 
open inner border. Positive torques arising from the inner 
disk are also stronger and To is reduced by 50%, i.e., the 
migration time scale is two times as long. 

When simulating a Jupiter-size body embedded in a 
disk with no initial gap, a density indentation is gradu- 
ally carved in. In order to skip the gap formation phase, 
an approximate analytical gap is sometimes imposed in 
the initial density distribution (see Kley 1999). We per- 
formed three computations adopting this choice. In these 
cases, a partial shrinking and refilling of the analytical 
gap is observed. Besides, material drains out of the inner 
radial border faster than it docs in our standard models 
(no initial gap). Hence, the inner disk depletion is in- 
tensified. With respect to standard models, we measure 
smaller accretion rates and longer migration time scales. 
Discrepancies in both quantities stay below 20%, after 
200 orbits. However, since the model outcomes indicate 
a tendency to converge as the evolution proceeds, a more 
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appropriate comparison should be made after a long-term 
evolution. 

4.5.1. Grid Resolution 

Hardly any hydrodynamic calculation is strictly reso- 
lution independent. Thus, for completeness we analyze 
in this section how our estimates on migration and accre- 
tion vary because of different hierarchy resolutions. Two 
tests are presented for each of the quantities Mp and 7d . 
They arc computed with the aid of grid systems G3 and 
G4 and then results are tested against those calculated 
with the less resolved hierarchy G2 (see Table 2). In 
this way we aim at checking finite resolution effects in 
the radial and azimuthal directions and, separately, those 
in the meridional direction. In fact, G3 and G2 have the 
same number of grid points in the vertical direction but 
Ai?[G3] = 0.82Aii[G2] and A<^[G3] = 0.75Av5[G2]. In 
the other test case (G4 against G2), R and ip gridding 
is unchanged while the number of latitude grid points is 
nearly doubled. 

With regard to the mass accretion rate (Fig. 9, left 
panel) difi^erences are below 20%. Though the number of 
grid cells in the accretion sphere is enlarged by a factor 
either 1.7 or 2, no systematic tendency seems to arise from 
this test. Something different happens to the total torque 
(right panel). In fact, while the increased resolution in 
the latitude direction does not play any considerable 
role, the larger number of grid points in R and ip causes a 
reduction of the total torque magnitude between 20 and 
30%. This is not at all unexpected. On one hand circum- 
stellax disk spirals are better captured by a finer gridding 
in the radial and azimuthal dimensions. Thus, Lindblad 
torques arc accounted for in a more accurate fashion. The 
Figure proves this to be especially true when the ratio q 
is small, because of the diminishing wave amplitudes. On 
the other hand, due to the vortical exponential drop of 
the density and the lack of temperature stratification, disk 
layers above the midplane do not contribute very much to 
To . A finer resolution along the vertical dimension cannot 
sensitively modify the total torque outcome. 

5. Discussion 

Here we devote some further comments to the differ- 
ences between accreting and non-accreting protoplanets 
and then to the effects of the vertical density structure on 
the gravitational torques acting on embedded objects. 

5.1. Pressure Effects in Protoplanetary En- 
velopes 

For the purpose of carrying out a local analysis in 
the vicinity of accreting and non-accreting protoplanets, 
we introduce a cylindrical coordinate reference system 
{O' ; I , , z} with its origin O' coinciding with the planet 
position and the z-axis perpendicular to the disk mid- 
plane. Hence, we will have that S'^ = l'^ + z'^. The lon- 



gitude angle i/> is counterclockwise increasing, and tp = n 
points toward the star. Supposing that the flow nearby the 
planet is stationary, neglecting fluid advection and viscos- 
ity, the Navier-Stokes equation for the radial momentum 
reads: 

«4_a$p i9p 
I ~ di ^ p dv ^ ' 

where is the azimuthal velocity component around 
the planet. Excluding the particular situation represented 
by a homogeneous sphere (eq. [6]), the first term on the 
right hand side of equation (17) is positive (see Fig. 1). 
Recalling equation (2), we see that the second term is 
proportional to the density gradient, which is negative, 
and therefore it reduces the centrifugal acceleration w'^jl. 
In Figure 10 we show some quantities, at « = 0, aver- 
aged over the angle %[>, regarding the same simulations ad- 
dressed in § 4.1.1 (Mp = 20 Me with <I>p = <i>l^). From 
the top left panel one can realize that the mean density is 
indeed higher in the non-accreting case (solid line) than 
it is in the accreting case (dashed line), as it weis argued 
in § 4.1.1 from the values listed in Table 4. In order to 
evaluate how much the pressure gradient affects the left 
hand side of equation (17) in both cases, we plot the av- 
erage of such quantity {{w^)/l) in the top right panel of 
Figure 10. The centrifugal acceleration is much smaller in 
the envelope of the non-accreting model (solid line) than 
it is in that of the accreting one. Such circumstance is 
a clear indication that the envelope is pressure supported 
in the first case. The behavior of the averaged velocities 
(ui^) and (wi) is shown in the two bottom panels. As 
expected in a pressure dominated fiow, the magnitude of 
both velocity components is smaller in the non-accreting 
model (dashed lines). 

5.2. Torque Overestimation in 2D Geometry 

Gravitational torques exerted by a three-dimensional 
disk onto a medium- or low-mass protoplanet are weaker 
than those generated by a two-dimensional disk. Miyoshi 
et al. (1999) state that the total torque 7d in 3D is 0.43 
times as small as that in 2D. Something similar was found 
by Tanaka, Takeuchi, & Ward (2002). Our fully non-linear 
calculations predict that low-mass protoplanets have a mi- 
gration rate, at least, an order of magnitude less in disks 
with thickness than they have in infinitcsimally thin disks. 
One of the main reasons for that relies upon the vertical 
decay of the density, as one can demonstrate easily with 
a simplified approach. 

Let tz be the ^-component of the gravitational torque 
exerted by a column of mass Tildldip, located at distance I 
from the planet (see § 5.1). The surface density is defined 
as S = J pdz. If /g is the force exerted by such mass 
distribution, projected on the equatorial plane, then we 
can write 

t, = RpfsSmi;. (18) 
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Fig. 9. — Left panel. Comparison of the accretion rate Mp as computed on grid systems with different 
resolutions. Squares indicate that the ratios Mp[G4]/Mp[G2] arc drawn, whereas triangles refer to the ratios 
Mp[G3]/Mp[G2] (see text). Right panel. The same type of comparison for the total torque 7d experienced by 
the protoplanet. Squares and triangles have the same meaning as before. At Mp = 10 M^, we compare results 
from models with the homogeneous sphere potential given in equation (6), when available (for this particular 
mass value, see the discussion in § 4.2). 



The ratio of tz [3D] to tz [2D] is therefore equal to 



/a [3D] _ 
/«[2D] E 



r 

J — C 



{P + «2)3/2 



dz. 



(19) 



Since /g[3D] and /9[2D] are coherent in sign, x is also 

equal to the ratio of |fz[3D]| to |fz[2D]|. In order to quan- 
tify this quantity, we can assume a Gaussian mass density 
profile with a scale-height H, which is appropriate as long 
as no deep gap has formed. Thus 



-f 



exp 



{P + ^2)3/2 



dz. 



(20) 



The ratio x function of I is plotted in Figure 11 
and it evidences how a two-dimensional geometry overes- 
timates the magnitude of gravitational torques acting on 
the protoplanet. In the limit P ^ 2^, x converges to 1, 
which proves that only torques arising from locations near 
to the planet (/ < H) are magnified. 

Though larger torque magnitudes do not necessarily 
imply faster migration speeds, they can favor a larger mis- 
match between negative and positive torques and there- 
fore shorter tm- 

6. Conclusions 

On the background of the numerical computations of 
disk-planet interaction presented in Paper I and II, in this 



paper we combine the full 3D geometry of a circumstellar 
disk with a nested-grid technique in order to investigate 
in detail fiow dynamics, orbital decay, and mass accre- 
tion of protoplanets in the mass range [1.5 M®, 1 Mt^]. 
Besides, we overcome the point-mass assumption by em- 
ploying analytic expressions of the gravitational potential 
derived from simple theoretical models of protoplanctary 
envelopes. Each of them applies to distinct physical situ- 
ations: when the envelope mass is negligible with respect 
to the core mass; when the envelope is homogeneous and 
much more massive than the core; when the envelope is 
fully radiative, and finally when it is fully convective. 

Through a series of 51 simulations, we inspect the 
evolution and differences of protoplanets represented by 
tlic aforementioned gravitational potentials. We analyze 
the behavior of both accreting and non-accreting objects. 
Furthermore, we evaluate physical and numerical effects 
due to our standard set-up of the models. The com- 
putations clearly show that to accurately determine the 
early physical evolution of planets three-dimensional ef- 
fects have to be taken into account. 
The main results of our studies can be summarized as: 

1. Above the disk midplane the flow is nearly laminar 
only far away from the planet. The region of in- 
fluence of the planet extends well outside the Hill 
sphere and its boundaries are marked by vertical 

shock fronts. Past the shock, matter is deflected up- 
ward and then downward. In some cases, a closed 
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Fig. 10. — Midplane quantities {z = ov = 7r/2) azimuthally averaged around the protoplanet for two 
simulations in which <I>p = and Af,, = 20 Af®. The sohd hnc belongs to a non-accreting model, the 
dashed line to an accreting one. Top-left: mass density. Top-right: centrifugal acceleration. Bottom-left: 
rotational velocity. Bottom-right: velocity component along the radial distance I (see § 5). All of the four 
quantities indicate that the envelope is mostly pressure supported in the non-accreting case, whereas it is mainly 
centrifugally supported in the other. 



recirculation is also observed. In the disk midplane, 

spiral waves around the planet arc not as strong and 
tight as they appear in two dimensions because of 
wave deflection in the vertical direction. 

2. In the mass range of their applicability, Steven- 
son's and Wuchterl's gravitational potentials pro- 



duce flow structures, close to the planet, similar to 
those determined by a smoothed point-mass poten- 
tial. Migration times and accretion rates arc alike. 
In contrast models with the (unrealistic) potential 
of a homogeneous sphere yield different dynamics 
though, as for tm and Mp, not much difference is 
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Fig. 11. — Ratio of the tree-dimensional to the two- 
dimensional torque exerted by a column of material 
lying at a midplane distance I from the protoplanet 
(see § 5). The upper x-axis is in units of the disk 
semi-thickness at the planet location: H = hRp. 

observed for Jupiter and Earth-size bodies. 

3. Since the numerical accretion procedure might be 
considered somewhat arbitrary, we ran several mod- 
els in which the protoplanet does not accrete at 
all. Non-accrcting models behave differently from 
accreting ones in a volume whose size is roughly 
comparable with the Hill sphere. Within this re- 
gion matter is pressure supported and thus a spher- 
ical envelope builds up. Except for the case Mp = 
20 M®, the total torque Td exerted by the disk is 
on the same order of magnitude as that measured 
in accreting models. 

4. According to Ward's theory (Ward 1997), the mi- 
gration speed settles to a constant value when the 
planet- to-star mass ratio g > 4x 10~*. Our numeri- 
cal results give a similar trend at a slightly different 
magnitude though. Most of these simulations pre- 
dict an inward migration except the one where a 
20 M0, non-accreting, protoplanet is involved. In 
the mass range [7 , 15 M®] migration speeds can 
be 30 times as slow as those predicted by Tanaka, 
Takcuchi, & Ward (2002) although, outside of this 
range, the agreement between our computational 
data and the type I migration by the same authors 
is remarkably good. We suspect that this surprising 
outcome may be caused by the onset of non-linear 
effects appearing around ten Earth's masses, which 
conspire to give such long migration time scales. If 



correct, this much slower inward motion may help 
to solve the problem of the too rapid drift of planets 
toward their host stars. 

5. In agreement with studies on planet formation (Bo- 
denheimer & Pollack 1986; Tajima & Nakagawa 
1997), the growth time scale shortens as the pro- 
toplanct's mass increases. The minimum is found 
at Alp = 20 Mq). Albeit the feeding process slows 
down as soon as angular momentum transferred 
by the planet to the surrounding material is large 
enough to dig a density gap. Then, at Mp w 1 Mi^, 
the accretion rate greatly reduces and the growth 
time scale becomes consistently very long. We 
present an analytical formula for the growth rate 
which may be useful for global studies in planet 
formation. 

6. As long as migration and mass accretion are consid- 
ered, two-dimensional computations still yields re- 
liable results when the mass ratio q > lO""* (Mp > 
30 M(g if Mi, — 1 AIq). In practice, 2D geometry is 
applicable whenever the Hill radius Rh exceeds the 
60% of the local pressure scale-height of the disk 
H. But for smaller masses three-dimensional calcu- 
lations have to be considered. 

The three-dimensional calculations presented here 
achieve a new level of accuracy by using a sophisticated 
nested-grid technique. This numerical feature allows a 
global and local resolution not obtained hitherto. How- 
ever, similar to all of the previous calculations, the models 
presented here have one principal limitation: the lack of 
an appropriate energy equation. Because of this, we could 
not couple the thermal and the hydrodynamical evolution 
of the system. If one wishes to do that in three dimen- 
sions, the energy equation has to include radiation and 
convective transfer. Yet, only with massive parallel com- 
putations one can hope to pursue this goal. Thereupon, 
the direction of future developments and improvements is 
already marked. 
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A. Interpolation Formulas for the Fine-coarse Grid Interaction 

The aim of this appendix is to furnish some algorithms useful for updating the values of scalars and momenta 
on a coarse grid with those computed on the hosted (finer) grid, in spherical polar coordinates. 

For the purpose, we indicate with the mass density to be interpolated on the coarse level. The interpolating 
values, on the finer subgrid level, are indicated simply as p{i,j,k). For the density average, the i-index varies 
between i and i + I and so do the other two indexes. 

Accordingly, U^, and are the coarse linear, meridional and azimuthal angular momentum. The finer 
values from which they are reset will be denoted as UR{i,j, k), U${i,j, k), and U,p{i,j, k), respectively. For the 
average of the radial vector component, the i-index varies between i — 1 and i + I, while the other two indexes 
vary between j (k) and j + I {k + 1). This is related to the locations where velocity components are defined 
on the mesh. In case of the meridional angular momentum, the j-index ranges from j — 1 to j + 1, while the 
others are i (fc) and i + 1 (fc + 1). For the azimuthal angular momentum, it is the fc-index to extend over the 
wider range. 

The finer grid coordinates are {Ri,9j,(pk), and the grid spacing is (Ai?, A^, At^) (see Fig. 12), which is 
constant. Since the grid has a staggered structure, scalars are volume-centered, i.e., p(i,j,k) lies at (Ri + 
AR/2, 0j + A0/2, (fk + A.(f/2), while resides at {Ri + AR, 0j + AO, ipk + Atp), because the linear resolution 
doubles from a grid to the hosted one. Instead, vector components are centered each on a different face of the 
volume element. For example, the radial component UB,{i,j, k) is located at (i?j, 9j + A9/2, (fk + Alp/2), whereas 
the coarse radial momentum is defined at {Ri, 9j + A9, (pk + A(p). The locations of the other components 
follow by similarity. 

The interpolation is basically a volume-weighted average. Eight volumes are necessary to carry out a scalar 
interpolation. Momentum interpolations require that twelve spherical sectors must be employed. Yet, since the 
metric in a spherical polar topology is independent of the azimuthal angle p, some of them actually coincide. 
In order to distinguish among the four volume sets, we introduce the notations V^p\ V^^\ V^^^ and V^'^\ 
according to the quantity to average {p'^, U^, Uq , and U^). These four sets differ because of the space metric 
and the staggered mesh. 

Once the correct elements have been identified, the coarse mass density can be replaced by 

In fact, sectors V^''^ are tp-independent, thus only four volumes enter this average. In a similar fashion, corrected 
momenta can be written, in a concise form, as 

where "E. = R, 9, and Lp. For computational purposes, a volume element is preferentially cast into the form 

V = (Ai?V3) (-Acos6') {A<p) . (A3) 
The four sectors V^^^ required in equation (Al) are the following 

V^''\i,j) = \ - Rl) [cos {9,) - cos {9,+,)] A<p 

yW(i + = 1 [(i?,+i + ARf - Rl+,] [cos {Oj) - cos (^,-+i)] A^ 

{i, j + = i (i?f+i - i?f ) [cos (^,.+i) - cos (%+i + A^)] A^ (A4) 

Vip){i + i,j + l) = ^[{R,+,+ARf-Rf^,] 

[cos {9j+i) - cos {9j+i + A0)] Aif. 
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Also for the radial and meridional directions, the denominator of equation (A2) reduces to 2 
though the summation includes six terms, this time. Therefore, the set of volume elements necessary for the 
interpolation of the radial momentum is: 



y(«)(z + i,j + i) 



[cos {9j) — cos {9j+i)] A(p 
[cos {9j) — cos (Oj+i)] A(fi 



??3 



- {Ri + Ri+i)^ 



[cos {6j) — cos (Oj+i)] Atp 



{Ri-i + Ri)^ — -Rf_i 



[cos {6j+i) - cos {Oj+i + AG)] Aip 

— [{Ri + Ri+i)^ — {Ri-i + Rif\ 
[cos (6*^+1) - cos [Oj+i + AO)] Aip 



Ri+i - - + Rt+i) 



[cos {9j+i) - cos (6*^+1 + AO)] A(p. 



(A5) 



The group of elements involved in the updating process of the meridional angular momentum Ug is: 



V^'\i + l,j-l) 



= ^ - Rl) {cos (^,_i) - cos [{Bj-i + ej)/2]} A^ 



3 (-Rf+i - -Rf) 

-cos[(^,-+^,-+i)/2]} 

i?^){cos[(a " ■ 



3 - 

{cos ^,0/2] . 

-e,+i)/2]-cos(0j+i)} A^ 



i [(i?,+i + ARf - 
{cos {0j-i) - COS + 0j)/2]} A>fi 
i [(i?.+i + AR)^ - Rl^] 

' -cos[(e, +e,+i)/2]} 



3 

{cos[{0j.,+ej)/2] 



^[{Ri+i+ARf-Rl^] 
{cos [{9j + ej+i)/2] - cos (^,+1)} A^. 



(A6) 



In order to perform the interpolation of the azimuthal angular momentum , the following set of volume 
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Ri-i 

















'k+l 



Fig. 12. — Sketch of a spherical sector composed of two coarse grid cells (thick lines) and the set of nested cells 
(thin lines). The discretized spherical coordinates (i?i, Oj, (fk) are relative to the fine grid, and so is the resolution 
{AR,A9,A(fi). Scalars are cell-centered on each grid level, whereas vector components are face-centered, each 
on a different face of the volume element. 



elements is required: 

V^^\z,j, k-1) = i {Rl, - i?3) [cos (9,) - cos (9,+^)] A^/2 

vM{i + l,j,k-l) = i [(i?i+i -h Ai?)3 - i?f+i] 

[cos(6lj) - cos(6'j+i)] Atp/2 
V('P){i,j + l,k-l) = ^{Rl,-R^) 

[cos {6j+i) - cos {9j+i + A9)] Ay/2 
VM{i + l,j + l,k-l) = ^[{Ri+,+ARf-R^^^] 

[cos {Oj+i) - cos {9j+i + A9)] Aip/2 
V^'^Hi,j,k) = 2W^)(i,j,fc-l) (A7) 
F('^)(i,j + l,fc) = 2V('^)(i,.7 + l,fc-l) 
y('^)(i + l,j,fc) = 2V^'^'>{i + l,j,k-l) 
V^'^\i + l,j + l,k) = 2V^'^\i + l,j + l,k-l) 

v^'^\i,j,k + i) = y(^)(i,i,fc-i) 

V'^'^\i,j + l,k + l) = V^'^\i,j + l,k-l) 
V^'<''>{i + l,j,k + l) = V^'^'>{i + l,j,k-l) 
V^'^\i + l,j + l,k + l) = y(^)(i + l,i + l,A;- 1). 

Equations A7 imply that the denominator of equation (A2) is also equivalent to 2 V^'^^i,j, k). 

Velocities arc retrieved from momenta and density. Fine-coarse interaction also involves the correction of 
momentum flux components across the boundaries between two neighboring grids. This is accomplished via 
time and surface-weighted means (for details, see Paper II). 
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